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SUMMARY 

In this paper vector potential and related methods, for the simulation of both 
inviscid and viscous flows over aerodynamic configurations, are briefly 
reviewed. The advantages and the disadvantages of several formulations are 
discussed and alternate strategies are recommended. The paper consists of the 
following sections 

Introduction 
Scalar Potential 
Modified Potential 

Alternate Formulations of Euler Equations 
Least-Squares Formulation 
Variational Principles 

Remarks on Iterative Techniques and Related Methods 

Viscous Flow Simulation 

Conclusions 


INTRODUCTION 

Most of the recent successful attempts to solve Euler equations are based on the 
primitive variable formulation where the governing equations are written in a 
divergence form representing the conservation laws of mass, momentum and energy. 
The main reason is the capability of correctly capturing flow discontinuities 
(shocks and wakes) with the help of artificial viscosity terms, explicitly added to 
the equations or implicitly built in the numerical schemes. Moreover, the 
conservation laws are basically a hyperbolic system of equations which can be 
easily integrated in time (using for example explicit schemes like Runge Kutta 
methods) to obtain the steady state solution. This fact is particularly attractive 
if unstructured grids are used for simulating flows over complex three-dimensional 
configurations. Such a calculation can be very efficient on present computers 
since it is fully vectorizable and indeed very impressive results were obtained 
using multigrid convergence acceleration techniques. 

Artificial viscosity may lead, however, to artificial vorticity, artificial 
boundary layers and separation. To minimize the contamination of the inviscid 
solution with such artificial viscous effects, higher order schemes are preferred 
and to avoid overshoots and undershoots near shocks, flux limiters are applied. 
Finite-difference schemes, with the total variation diminishing property, are 
proposed. Because of their diagonal dominance, relaxation methods are applicable and 
hence they are used as smoothers for multigrid techniques. The limiters, however, 
are highly nonlinear and some difficulties are expected, such as nonexistence and 
nonuniqueness of the discrete solution of the conservation laws, as well as slow 
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convergence and limit cycles. Research efforts will continue to construct high 
resolution schemes particularly to capture contact discontinuities, and successful 
applications of multigrid acceleration techniques will follow. 

The purpose of the present paper is to examine alternate approaches, particularly 
for aerodynamic applications. The rationale, here, is very simple; far away from 
the body, the disturbances are small and the flow can be adequately described by a 
single equation, with a single unknown (the potential function) and a single 
parameter (Mach number); the potential formulation is valid, as long as the 
vorticity in the field is negligible. 

One obvious choice is the zonal approach, where the potential formulation is 
restricted to the irrotational flow zone while the Euler primitive variable 
formulation is used for the inviscid rotational flow zone. One problem with this 
approach is the matching of the two local solutions along the interface or in the 
overlap domain of the two zones. The quality of the solution and the convergence 
of the overall calculations will depend on the numerical treatment of the 
interface problem. 

Another approach is to augment the potential formulation with the rotational 
effects. This can be done using the Helmholtz decomposition theorem, where the 
velocity vector can be represented as the sum of a gradient of a scalar potential 
function and the curl of a divergence free vector. The correction to the potential 
formulation automatically vanishes, when the flow is irrotational. The problem 
with this formulation is the implementation of proper boundary conditions for the 
corrections particularly for multiply connected domains. 

Related to this approach is the use of a variational principle and the Clebsch 
transformation. The connections to least-squares formulations are also delineated. 

Extensions of these methods for viscous flow simulations are discussed and the 
relation to some viscous/inviscid interaction procedures are depicted. 

In the following, recent work on scalar potential methods are reviewed first and 
then details of the above approaches are studied. Finally, some concluding remarks 
are drawn. 


SCALAR POTENTIAL 


In the recent literature, one can find excellent reviews of methods of solution 
of the potential equation (See for example Caughey (1) , Holst (2) and South (3)). 
In this section, only the most recent developments will be reviewed. 


Steady inviscid flows, with constant entropy and total enthalpy, can be 
described by a potential function satisfying the continuity equation in 
conservation form 

V.(pV^) = o (1) 

where 

1 

e - (l-V M* ( |v*| -1)) 1_1 (1') 

and M is the free stream Mach number. 

00 
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The full potential equation (1) is nonlinear and of mixed type (elliptic in the 
subsonic region and hyperbolic in the supersonic region). It admits expansion as 
well as compression shocks. Honest discretization of equation (1) leads to a 
discrete system of equations with the same properties. One way to exclude 
expansion shocks is to use artificial viscosity methods. No vorticity can be 
generated due to the artificial viscosity in such calculations since the velocity 
field is calculated as a gradient of the potential function. Moreover, for pure 
subsonic flows, no artificial viscosity is needed and the potential solution 
should satisfy the Euler equations. On the other hand, to obtain the solution 
of Euler equations, artificial viscosity is usually needed even for pure subsonic 
f 1 ows . 

The artificial viscosity can be introduced in the numerical scheme in many ways. 
Osher et al . (4) introduced a flux biasing scheme which, unlike the retarded 
density scheme, does not allow for expansion shocks in the discrete solution of 
equation (1). This scheme is first-order accurate in the supersonic region, while 
the artificial viscosity is switched off in the subsonic region. Shankar et al. 
(5) used it successfully for time -dependent calculations. Volpe and Jameson (6) 
applied a second-order version of the same scheme in their multigrid code and 
obtained impressive results. 

Recently, Mostrel (7) introduced a second-order accurate scheme, with no 
subsonic/supersonic switching, and proved global linear stability, total variation 
diminishing (with flux limiters) and discrete entropy inequality. He also 
introduced a time splitting algorithm for the 2D low-frequency transonic small 
disturbance equation. In this area, one should mention the unpublished work of 
Catherall* who introduced an approximate factorization scheme, AF2, which consists 
of two factors only, even for three-dimensional flows. 

Various methods are used to accelerate the convergence of potential flow 
calculations. Besides multigrid, generalized conjugate gradients, generalized 
minimal residuals (GMRES) and extrapolation procedures can be applied on computers 
with large memories. For two dimensional (and axisymmetric) problems, direct 
solvers based on banded Gaussian elimination or nested dissection are nowadays 
feasible and quadratic convergence can be obtained using Newton's method (8-9). 
For three dimensional problems, block relaxation and domain decomposition 
procedures are needed. Some of these algorithms are suitable for parallel 
processors. 

Few years ago, many three-dimensional codes were developed for potential flow 
calculations. There are two common mistakes in most of these codes. The first is 
to use a two-dimensional vortex as a far-field boundary condition for each spanwise 
station. The second is to assume that the flow leaves the airfoil section in a 

direction bisecting the trail ing-edge angle as in two-dimensional problems; this 

is not true for general three-dimensional lifting flows. In all these codes, the 
wakes are fixed and are not allowed to adjust with the flow. No effort has been 

spent to correct these codes and it seems there is no interest to do so! In 

Europe, a finite-element code based on optimal control formulation was developed at 
the same time to simulate flow over a complete aircraft. It is not clear, however, 
how the vortex sheets are treated in such calculations. 

Recently, an intensive effort was launched at Boeing to develop a similar code 
(TRANAIR). Some of its innovative ideas are discussed in ref. (10). (An 
interesting test cas e is a wing in an incompressible flow where panel methods can 

*Catherall , D., private communication, July 1985. 
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predict correctly the induced drag. Moreover, theoretical results are available 
for certain configurations.) 

Finally, no progress has been made to further understand the nonuniqueness 
problem (ll) of the lifting potential flows in the transonic regime. Since no 
multiple solutions of the Euler equations for the lifting airfoil in the transonic 
regime have been reported, one tends to assume that the potential model is 
inadequate. It seems, however, that acceptable results can be obtained if boundary 
layer interactions are taken into account. After all, the potential flow 
approximation is only meaningful outside the viscous layers. It seems also that 
the viscous transonic equation (with uniform viscosity) has a unique solution but 
no proof is available. 


MODIFIED POTENTIAL 

In general, due to the absence of entropy and vorticity effects, shocks in 
potential solutions are stronger. In this section, approximate methods to account 
for these effects are discussed. 

It is argued that the vorticity effects are higher order. In fact, for straight 
shocks there is an entropy jump, but there is no entropy gradient downstream of 
the shock and, therefore, no vorticity is generated (it is well known that 
vorticity is related to the curvature of the shock). If the vorticity is 
neglected, the flow downstream of a shock can be also described by a potential 
function. The flow, however, is no longer isentropic. In this model, the entropy 
jump is accounted for by modifying the shock point operator. It turns out that the 
equation downstream of the shock is unaltered (since entropy is constant along 
streamlines). This approximation is completely equivalent to fit a Rankine- 
Huognoit shock in potential flows. Special treatment of wakes may be also 
required (12), (13). 

The next step is to construct a simple approximation for the vorticity effects. 
This can be achieved by augmenting the velocity due to the potential field by a 
rotational increment due to the entropy gradient. In two dimensions, an ordinary 
differential equation is solved for the rotational component, u, namely 


«... d( aS/ R) 
an d* 


( 2 ) 


For more details see ref. (12). 

An approximate solution of equation (2) is simply 



pq 


aS /R 

7M* 

00 


u 


( 2 ') 



Although the rotational correction is formally higher order than the correction 
due to the modification of the jump condition across the shock*, it is not 
recommended to neglect u, particularly for lifting airfoils. 

The reason is clear from examining the work of Klopfer and Nixon (14). Across the 
wake, the static pressure must be continuous while the total pressure, in general, 
is not. Therefore, the tangential velocity component must jump. If the flow is 
presented only by a potential function, the jump of the potential across the wake 
will grow linearly with the distance from the trailing edge (i.e. in the far field, 
r-«! ) . 

On the other hand, this nonuniformity does not occur if the potential field is 
augmented by the rotational component. In this case, both the entropy and u jump 
across the wake while 4> does not. This can be shown by expanding the static 

pressure formula (assuming the wake is aligned with the x - axis) 


P 


(1 



M„ 2 ((* x + U ) 2 




(3) 


The contribution of the entropy term cancels the contribution of the rotational 
component and therefore <t> must be continuous to guarantee continuous static 
pressure across the wake. 

For the three-dimensional flows, the calculations of the rotational velocity 
components are more complicated as will be discussed in the next sections. As a 
crude approximation, the two-dimensional formula can be used in each spanwise 
station. Dang and Chen (15) used instead, the following formula 


A** ft (4) 

2 

-rM (n cosa + n sina) 

® A Jr 


where ^ “ n x ^ x + n v normal to the shock surface. In our opinion, the extra 
computational erfortr irr the calculation of ft is not justifiable and the assumption 
of a locally normal shock is consistent with the other approximations made in the 
derivation of the rotational velocity component formulas ((2') as well as (4)). 


p 

* From Prandtl relation, (1+u, ) - (l+u«) = a* , the correction to the shock jump 
condition to allow for entnapy generation is second order, while u is third 
order. 
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ALTERNATE FORMULATIONS OF EULER EQUATIONS 


Two-Dimensional Flows 

In this section, no small disturbance approximations are assumed and the 
treatment of the exact inviscid equations is discussed. It is instructive to 

consider first a simple case of two-dimensional inviscid incompressible steady 
flow. The conservation laws are usually written in terms of the continuity and the 
momentum equations. On the other hand, the following equations are completely 
equivalent: 


u 


x 


+ v. 


0 


' u y + v x “ w (*) 


* + 2 (u 2 + v 2 ) * H (*) (5) 

where u> and H are constants along streamlines. Equation (5) represents 
conservation of mass, vorticity and total enthalpy. The first two equations are 
linear in u and v; the pressure is decoupled and is obtained from Bernoulli's law. 
If equation (5) is used , contact discontinuities (wakes) must be fitted and no 
artificial viscosity is needed in such calculations. 

The exact equation for two-dimensional inviscid compressible flows is 
w t + F x + G y = o (6) 


where W * ( p , pu, pv, pc)^ and F and G are the corresponding flux vectors. The 
Jacobian of F has the eigenvalues u + a, o, o, and similarly the Jacobian of G has 
the eigenvalues v + a, o, o, where a is the speed of sound. Four modes are 
identifiable: two acoustic modes, entropy mode and vorticity mode. In fact, for 
steady flows, the above equations can be rewritten in the following nonconservative 
form: 

(a 2 - u 2 ) u x - uvv x - uvUy + (a 2 - v 2 )v y « o (6a) 

-u + v x - « (6b) 


DaS/R _ 
Dt “ 


(6c) 

(6d) 


where 


Dt 


is the substantial derivative. 


The vorticity u> is calculated from the Crocco's relation 

D dASZE . , dH 
P d* p dv 


G 0 


(7) 
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The speed of sound is given by 


where 


2 

a 


(7-l)(H 


2 



( 8 ) 


2 2 
u +v , 


and the density and the pressure relations are 


/M 2 a " AS y n 

P * (M* a ) * e / R 


( 9 ) 


P ■ (" 2 a 2 ) 1 ' 1 • e ' iS /R /7H 2 (9') 

aS/R must have the proper jump across the shock to satisfy the Rankine Hugoniot 
relations. 

It is clear from equations (6a) and (6b) that the corresponding characteristics 
have the same form for rotational and irrotational flows. 

A shock fitting procedure must be used with this nonconservative formulation. In 
refs. (12), (16), the author used the Prandtl relations across the shock 


n „ _ a *2 ail „2 

q ln q 2n " a ‘ 7 +l q t 


(10a) 


q tj " q t 2 * q t ( 10b) 


For flows with constant entropy and total enthalpy, the vorticity vanishes and 
the irrotational motion can be represented by a potential field. The weak solution 
admitted by the potential equation (1) reflects conservation of mass under the 
isentropic assumption. On the other hand, in the modified potential formulation, 
the Prandtl relations are forced across the shock, which implies that the mass, 
including entropy effects, is conserved. To fully account for the entropy effects, 
a rotational velocity field due to the vorticity generated by the shock curvature 
must be added to the potential field; therefore, it is natural to use the 
decomposition 

u = <f x + u (11a) 

v * <t> y + v (lib) 

Obviously, such a decomposition is not unique; the rotational velocity components 
are not independent. A feasible constraint is 

u x + v y - o (12) 
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Such a constraint is automatically satisfied if a perturbutive stream function * is 
introduced, such that 


u - ♦ and v - - * (13-14) 

J * 

Thus, the governing equations for 4 and * are 

(p* x ) x + (p* y ) y - - (p* y ) x + (p* x ) y (15a) 

*xx + *yy w (15b) 

Equations (12) & (13) are solved in ref. (13). Chaderjian and Steger (17) used a 
similar approach for the lifting airfoil problem, where both 4 and * are chosen to 
be continuous across the wake. 


Another choice is 


u = * x + *v /p 

(16a) 

v = V V' 

(16b) 

In this case the governing equations are 


(p^ x ) x + (*^ y ) y = 0 

(17a) 

( Wy + (i x/p>y ‘ ‘ “ 

(17b) 

With proper boundary conditions for *, ^ can vanish 

identically. The 


In 


ref. (19), Papailiou et al . used the decomposition 


pu * <f> x + 5 y 

(18a) 

pv . * y - * x 

(18b) 

The corresponding equations are 


*xx * *yy ' 0 

(19a) 

and (V,)y ♦ (V,»x ■ ' " 

(19b) 


In all these cases, a partial differential equation for the correction must be 
solved. Other decompositions require _only the solution of ordinary differential 
equations. For example, the choice of v * 0 in (lib) leads to 

(p* x ) x + (p* y ) y = - (pu) x (20a) 


and u y • - u (20b) 

For certain grids, aligned locally with the flow, this decomposition is 
obviously useful. Another example is based on a multiplicative correction, 

u « A* x (21a) 

v - H y (21b) 
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The corresponding equations are 

(p* x ) x + (/^y )y * 0 (22a) 

*x Wx = ' w (22b) 

where p * Ap (23) 

The general Clebsch representation has both additive and multiplication 
corrections, as will be discussed later. Based on this discussion, the full stream 
function formulation is the most recommended one. It is not widely used for 
transonic flows, however, because of the difficulties associated with the nonunique 
relationship between the density and the flux. The same remark holds for 

axi symmetric flows (including swirl). 


Three-Dimensional Flows 
The conservation laws are given by 

w t + F x + G y + H z = 0 ( 5 ') 

The eigenvalues of the Jacobians of F, G, and H are u ± a,o,o,o, v ± a,o,o,o and 
w ± a, o,o,o. Accoustic, entropy and vorticity modes are still identifiable but 
the problem is more complicated since the vorticity has in general three nonzero 
components. 

For steady flows, the governing equations can be rewritten in the form: 

v . pq = o (24) 

vxq*S (25) 

Assuming p and Z are given, equations (24) and (25) are four equations in three 
unknowns. There is, however, a relation between the vorticity components 

v • u ■ o (26) 

The Crocco relation is used to calculate two of the vorticity components 

q x Z = - TvS + vH (27) 

The streamwise vorticity component can be calculated ^ using equation (26). 
Alternatively, taking the cross product of (27) with q, one can obtain the 
following formula for the vorticity (see refs. (21), (22), (24)) 

Z = A pq + x (Tvs - VH) (28) 

Equation (26) leads to an equation for \ 
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pQ* VA 


x (TvS - vH) 


(29) 



The entropy and the total enthalpy transport equations are 

q*vS = o (30) 

q*vH = o (31) 

While H and x are continuous across a shock, S jumps. The speed of sound, the 
density and the pressure are obtained from equations (8), (9) and (9 1 ), where q 2 = 

Idl- 

Using the Helmholtz theorem, <t can be decomposed into the gradient of a scalar 
function plus the curl of another vector 

q = v<f> + vx* (32) 

The first term on the right hand side of equation (27) _is curl free, while the 
second term is divergence free. A feasible constraint on 5 is 

v • ♦ * o (33) 

Hence, equations (24) and (25) reduce to 

v . (pv<t>) = - v.(p vxit) (34) 

v 2 * = - Z (35) 

The correction to the potential solution becomes a major effort, it requires the 
solution of three Poisson's equations in three dimensions. 


The boundary condition for the potential problem is 

n*v^ * - n*vx* (36) 

where n is normal to the solid surface. Two linear combinations of the stream 
functions can be kept constant on the boundary surface, and equations (33) can be 
used to solve for a third linear combination. 

For example, in orthogonal coordinates, two components ^of can be chosen to be 
constant in the plane tangent to the body, hence n*vxl = o, while the normal 
derivative of the third component vanishes. The boundary conditions in generalized 
curvel inear coordinates are given in refs. (20), (21), for simply connected 
domains. 


In a series of interesting papers (22), (23), (24), Dabaghi and Pironneau studied 
the vector potential method and obtained two- and three-dimensional finite-element 
solutions for transonic flow problems. The following variational problems are 
considered. For a given vector field q, find <t> and v such that: 
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(37) 


(V*,vw) - (q,vw) 

(vx», vxv) + (v*l,v.v) 

- (q,vxv) (38) 

where w and v are proper weighting functions (vxn| * o, J v»n d 7 - o ; r is the 
boundary of a simply connected domain). It can 1 be shown rigorously that the 
Helmholtz decomposition is unique 

q = v<£ + vx* and * o (39) 

The variational form of the vector identity 

vxvxi = - A? + (40) 

is given by 


(v*, vv) = (vx*, vxv) + (v.|, v* v) + 2 J r ^ d 7 (41) 

where R is the mean radius of curvature of r. Therefore, equation (38) for * is 
equivalent to 


- At - w 


M 

an 



o on r 


*x n| r = o (42) 

In general, the three components of ♦ are coupled through the boundary 
conditions. The formulation has been extended to multiply connected domains by 
Dominguez (25). Across the wakes, the potential function jumps and the 
circulations (the difference in potential values) remain constant, while the normal 
component of the gradient of the potential is continuous. If the wake surface is 
denoted by s. and the circulation by a. then 

J J 



The following constraints must be imposed on 

f tf*n d 7 = fi. (44) 

r i 

where n- are constants. The decomposition is no longer unique. It is unique, 
however, when n . are given and a. are adjusted so that 

* J 

v<f> + vx* = q on z (45) 

The extra difficulty for nonsimply connected domains stems from the fact that 
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0 


(46a) 


vxq * o , v«q * o , q* n | - 

r 

does not imply q * o (see ref. (26)). 

Unlike the primitive variable formulation of Euler equations, a special 
treatment of the Kutta condition (pressure is continuous at trailing edges) is 
required. 

Another choice for the velocity decomposition is 

q « v* + (vx*)/p (46b) 

With proper nonhomogeneous boundary conditions for *, <f> can vanish identically. 
This case is simply an extension of the stream function formulation to three 
dimensional flows. The equation for * is 

vx(vx*)/p - Z (47) 

and the boundary condition is 

n*vx*| r - q*n (48) 

Hirasaki and Heliums ( 27 )-( 29 ) and Richardson and Cornish ( 29 ) studied a similar 
problem. A partial differential equation on the surface of the boundary must be 

solved. The corresponding variational formulation is given by Dabaghi and 

Pironneau (24). The problem is to find a function g such that 

5 x n| r = g (49) 


implies equation (48) . 

It is shown that g = vf where f| is the unique solution of a Laplace-Beltrami 
equation 


s ( 

r 


2£ M 
asj aSj 


+ 


d$2 3S 2 


) d 7 


- f q*n wd 7 

r 


(50) 


where w is a weighting function andU^, S 2 ) is a set of an orthogonal local 
coordinate system on r. 

In ref. (20), Papal 1 i ou et al . used the decomposition 

P q = v<t> + vxl (51) 

and v«* = 0 . 

Therefore, the equations for <f> and * are 

v 2 ^ = 0 (52a) 

v 2 * = - v ^ x (^d) + pw (52b) 

p 
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Lagging the density in the right hand side of equation (52b) may lead to 
convergence difficulties for transonic flow calculations. 

In all these variations, the stream vector * has three components. In ref. (30), 
another method is proposed. One of the components of * is chosen to vanish 
identically, and viscous flows over a three-dimensional trough are calculated. The 
application to inviscid transonic flow calculations is given in ref. (31). 

Giese (32) proposed also a representation of the flux vector in three dimensions 
using two scalar stream functions (« and 6) 

pq « v*xv0 (53) 

The continuity equation is automatically satisfied since v*pq = v*v*xv*=o. In this 
formulation, the body is a stream surface. The equations for * and e are obtained 
from 


vx(v«r xve/p) = Z (54) 

No numerical solutions are reported in the literature for transonic flow problems. 
In ref. (33) an application to hypersonic flows around a body of revolution at an 
angle of attack is discussed. 

Recently, Rose (34) proposed to solve the Cauchy-Rieman equations using a finite 
element scheme where a three-dimensional potential solution in the element and only 
two dimensional stream function solutions on the boundary faces of the element, are 
required. 


LEAST-SQUARES FORMULATIONS 

In ref. (35), Fasel proposed to solve the incompressible flow equations using a 
ve^ocity/vorticity formulation, where the pressure is eliminated. The equations 
v.q = o and vxq = w are replaced by 

v 2 q = - v x Z (55) 

For viscous flow calculations, q = o on the solid surface and three Poisson's 
equations for three velocity components are solved, while Z is obtained from the 
vorticity transport equation. (Recently, Rose et al. (36) and Osswald et al. (37) 
solved the system of the first order equations directly). 

It seems that, unlike the stream function formulations, there is no difficulty 
with the boundary conditions. This is not true in general, since conservation of 
mass is not explicitly imposed and careful treatment of boundary conditions is 
required, particularly for inviscid flows. 

In ref. (38), the author proposed a least squares formulation with a systematic 
treatment of the boundary conditions. For the continuous problem, the following 
functional is minimized 

I (q) = Ja(v.q) 2 + |vxq-w| 2 dn (56) 
n 

where a is a Lagrange multiplier. The functional I is discretized first on 
structured (or unstructured) grids. Minimization of the discrete version of I, 
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with respect to the unknowns at each node leads to a conditioned system of 
algebraic equations. 

To demonstrate the relation between the least- squares and the vector potential 
formulations, we consider first the two-dimensional incompressible flow case. The 
equations, in general, are 


u + v 
x y 

-u + v 
y x 


= s 


<j> 


(57) 

(58) 


where the s term in equations (57) represents sources or sinks in the field. 
Equations (57) and (58) are written in the form 

L 0 * 0 < 59 > 

The least-squares formulation leads to the second order equations 

1*1 (J) - L*(*) (60) 

* 

where L is the adjoint operator. 

Alternatively, new variable <t> and * can be defined by the equation 

1*0 ■ 0 ( 61 ) 

Hence, equation (59) becomes 

(62) 


LI* 6 ■ (*> 


0) 


The variables <f> and # represents a potential and a stream function. Notice the 
right-hand side in equation (60) is differentiated while in equation (62), it is 
unaltered. 

The extension to compressible three dimensional flows is given in ref. (39). In 
cartesian coordinates, the equations are 

(pu) x + ( P v) y + = s 


V u y 
* w x + u z 


w. 


= w-2 


= (4)« 


= w. 


Equation (63) is written in the form 


a. 


v 

3 X 

0 


-a. 


V 

0 

■ 3 



' u ' 


s 


V 

- 

w 3 


w 


w 2 


W 1 




(63) 


(64) 
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The variables * *2 an( * *3 are introduced via the adjoint operator. 


■-P8 X 3y -> z 0- 

-p a„ -a„ 0 a. 


1 

CO 

1 

m 

’ pU ' 
PV 

y x z 

y d z 0 5 x " 3 y 


*2 

.*1. 


pu> 


or, - *>q = pV<f> + vxtf (66) 

If s = 0, <f> can identically vanish. Furthermore, if one of the ♦ components is 
chosen to be zero, a two-component stream function formulation can be obtained. 
Therefore, vector potential formulations are special cases of least squares. 


VARIATIONAL PRINCIPLES 


Unsteady Flows 


When Lagrangian coordinates are used, the equations of motion can be derived from 
Hamilton's principle (that the difference of kinetic energy and potential energy be 
stationary). In the Eulerian description, Seliger and Whitham (40) considered the 
functional 


rl <2 '«? ' '•» + * <!t + 


'<[ 


.ausx + ? lpq i S) ) 

+ * \ at ' 


* f ^ dx dt 


(67) 


where <f>, r/ and 0 are Lagrange multipliers. 

The variations with respect to q^, p, S, a lead to 


M. 


2jl 


fiq i : q, * “ + S “ + a 


ax i 


ax. 


SA 

ax i 


Sp : h 


_ 1 Q 2 _ fii _ c Du 
" 2 q i Dt Dt 5 Dt 


5S: 


Dn 

Dt 


- T 
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( 68 ) 


6 a 


. fifi . 

• Dt 


o 


where 


e + 


x_ E 

7-1 P 


The variations of and p reproduce the constraints 


a. a ( PQj ) 

s * : ft + fiTj 0 




( 68 ') 


The coordinate a does not change along a particle path and the terms avp in the 
velocity allow the representation of an initial vorticity to be separated from 
those produced by subsequent entropy gradients. In this representation the 
vorticity is given by 


w * VS X Vr; + Va X Vj0 (69) 

Seliger and Whitham (40) derived also a simplified form of the variational 
principle (67). Upon integration by parts the integrand L in (67) is equivalent to 

I _ i n 2 a c fin Cfi nn\ 

L - 2 p q - pe - p J) t - pS Dt - pa Dt (70) 

(all variations are taken to vanish on the boundary R) 

From Equation (68), 

L * ph - pe * p (71) 

Therefore, one can start directly with the functional 

JT P dX dt (72) 

R 

where p * p (h,S) 

Ph = p* Ps * ‘ pT 

Introducing a Clebsch representation for the velocity 

q - v<f> + SVfj + avp (73) 
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where 


(74) 


h 


- <j>. - Sq. - a p. 




the variations of the integral of the pressure with respect to <f > , S, fj, a give 
exactly the corresponding equations of (68) and (68'). 

While the form in (72) is far from Hamilton's, it can be considered as a 
generalization of the well known Bateman principle of isentropic flows. It is 
noticed, however, that a, p are not uniquely defined by (73) since any perfect 
differential may be added to <f> with consequent changes in a and p. From the 
equations of motion, it can be shown that 

ft (n, e} . 0 
a (x, y, z) 0 


1 2 

where n * h + £ q. + + S»? t + ap^ (75) 


Therefore, n * n (a, p , t) and a and p can be obtained from the Hamiltonian form 


Da £0 Qfi £0 

Dt * " ap ’ Dt " da 


(76) 


It is sometimes advantageous to retain n * o and to use (75) in place of (74). 
The variational principle (72) is not modified but the variations with respect to 
q and p then give (76) as required. 

In fact, some of these ideas are very old. The velocity representation in (68) 
was first introduced by Glebsch in 1859, for the case of incompressible flow. 
Lamb (1932) considered a baratropic flow (p = p(p)), where any velocity can be 
represented by 

q = + o7 p 


and the momentum equations can be integrated to yield 


where 


S f * 2 q i ' • *X - °^t 
and 


Da 

Dt " 0 


M m 0 

Dt 0 


The alternative Hamiltonian form of a and p equations was first introduced by 
Stuart in 1900. 

Recently, Buneman (41) showed that the flow equations can be put in a Hamiltonian 
form, with benefit for numerical schemes. He considered the isentropic flow case 
with Clebsch variables 
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pq * - a V/i - pv<f> 


( 77 ) 


Thus, the vortical and irrotational parts of the flow are represented in a 
symmetric manner. Note that there would be a similar third term if entropy changes 
are taken into account. In terms of p , 4>, a and p, the following functional 
derivatives are obtained 


as 

dip 


v»pq 


at 


m 

dp 


c ip 1 ' 1 - 2 ( J ) 2 ( v ^) 2 + 2 


M 

at 


SSL _ 

da ' 


q.7 M 


sbk 

at 


dp. 


v.aq 


a° 

at 


where 


n 


C P 7 + 


(aVp + pV<j>) / 


(78) 


and c is a constant. 

The first equation is conservation of mass, and the second is a generalized 
Bernoulli equation. The third and the fourth equations are statements that the 
vortex parameters p and a/p follow the flow. With this Hamiltonian form, one can 
stagger data for densities and potentials in space as well as in time. Also, 
densities can be updated by leap-frogging over potentials and vice-versa. No 
numerical solutions based on this formulation have been reported. 

On the other hand, Ecer and Akay (42) were successful to calculate rotational and 
transonic flows using a variational principle similar to (67). They considered 
the functional 


If [(| eq? - pe) + * { + 


,K) 


ax i 


+ 


Da. 


p Dt 




dX dt 


+ rj f dr dt + SS f fi dr dt 

trj 9 tr 2 p k 


(79) 


where rj and ^ are those portions of the boundaries where the following flux type 
quantities are specified. 
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f 4 m ^* n 


f O k ' '“k 

The variations of (79) with respect to a and p give 

Da. 

sp \C p dT = 0 


, . !!^k _ dH_ . E 

Sa k‘ p Dt ■ p da,, R da. 


(80) 

(81) 


(82) 


Another choice is also suggested in ref. (42). The total enthalpy is employed as 
a primary variable to replace the material coordinate a in (67). The proposed 
variational principle is, however, incorrect, since the energy equation for 
unsteady flow is 


DH _ 1 &> 
Dt ' P dt 


(83) 


and H does change along a particle path! 

Steady Flows 

Roberts (43) used directly the following decomposition for steady flows 

q = % + v* + (S-S Q ) Vf? + (H-H Q ) vp (84) 

where the equations for r? and p are given by 


Bn _ 

Dt " ' 


T 


M m j 
Dt 1 


(85) 


For multiple shocks, multiple fields must be introduced and the contribution 
to the velocity is given by 


k 

E (S t - Vr;^ (86) 


Similarly, multiple p fields must be introduced for multiple wakes. On the other 
hand, for isoenergetic flows, p vanishes identically. Grossman (44) calculated 
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successfully supersonic conical flows, with strong rotationality due to shocks, 
using such a decomposition. 


For a thin body in a uniform flow, using small disturbance approximations, 
equation (85) becomes 

“ M - - T <«'> 

substituting equation (85') into (84) yields 

Q - 1 + v<f> AS/ R ^ i 

« 1 + v* - && i (84') 

7 M c 

* oo 

The rotational component in equation (84') is identical with equation (2') which 
is derived based on Helmholtz decomposition. 


A modified form of Clebsch representation was introduced by Hirsch et al. (45). 


q - v<f> + vS + #2 


(87) 


Substituting equation (87) into the Crocco relation, one obtains for arbitrary and 
independent entropy and rothalpy gradients, the two equations for and *2 


It = 


T 


D*2 

“Dt 


- 1 


( 88 ) 


A simplified representation can be obtained if a unique relation between S and H 
exists in the inlet flow field. In ref. (46), the authors discussed applications 
to rotational internal subsonic flows in ducts. 

It is clear from the above discussion that the computational effort required for 
the variational principle formulation is less than that needed for the 
implementation of least squares or vector potentials. Thanks to Clebsch, we have a 
generalized Bernoulli equation for rotational, nonisentropic steady/and unsteady 
flows . 


REMARKS ON ITERATIVE TECHNIQUES 
AND RELATED METHODS 

A common problem in all the various formulations discussed so far i s the 
solution of a nonhomogeneous potential equation, where the right-hand side 
represents the rotational effects. A crucial point in developing an iterative 
technique to solve this nonlinear problem is to account for the mixed type nature 
of this equation. 
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Moreover, the characteristics of this equation should be based on the total 
velocity and not on the irrotational part only. 

Techniques based on lagging the density as well as the rotational components are 
not reliable for transonic flow calculations. For example, the potential equation 
for two dimensional flows may be written in the form 

(A 1 * x ) x + (A 1 * y ) y - - (A 1 u) x - (A'v) y (89) 

where n is the iteration count. Obviously, one should expect convergence 
difficulties since an asymmetric operator is replaced by a symmetric one. In 
practice, a crude solution may be obtained, if excessive artificial viscosity is 
used. To construct a scheme, independent of the artificial viscosity, equation 
(89) may be replaced by 

^2 ( (a 2 - u? ) s * xx - 2 uv 6 * xy + * v2 )** yy ) - - (/>u) x ' (*> v ) y ( 9 °) 


It is not important to keep the left-hand side of equation (90) in conservative 
form and there are well-known methods to solve such a mixed type equation. 
Note, the potential correction will not affect the correction to the vorticity 
equation and ideally, the rotational components can be constructed such that they 
conserve mass. In fact, the Helmholtz decomposition can be applied to the 
correction rather than the original variables. To demonstrate this application of 
the Helmholtz theorem, consider the simple Cauchy-Riemann problem 


u + v 
x y 


u. 


- to 


-y x 

One can construct the following two-step iterative technique 


(91) 


step 1: 


step 2: 


‘ u x + iv y " ' < u x + v y * s > 

5U y - 5V X = 0 
*u x + 6V y = 0 

fiu y * 5v x " * ( u y * v x + 


(92) 

(93) 


The correction in the first step is curl free, and according to the Stokes 
theorem, it can be represented by a potential. On the other hand, the correction 
in the second step is divergence free and according to the Gauss theorem, it can 
be represented by a stream function. (For a general three-dimensional problem, 
the correction in the second step can be represented as a curl of a divergence 
free vector.) 


One can use Helmholtz decomposition directly on the discrete velocity field. A 
similar idea is behind the distributive relaxation discussed by Brandt in ref. 
(47). In such a technique, each discrete equation is satisfied in its turn by 
distributing changes to the several unknowns appearing in the equation in a 
specific manner; the main property is that in relaxing one equation, all the 
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residuals of the other equations are kept unchanged. Application to the Cauchy- 
Riemann problem on a staggered grid is discussed in ref. (47) and it is completely 
equivalent to a discrete Helmholtz decomposition. 

Brandt also considered Stokes problem, where the above strategy is slightly 
modified. The equations are 


v.q « s 

(94) 

-vq + vp « f 

(95) 


First, lagging the pressure, the velocity components can be updated from the 
momentum equations. In the second step, the velocity must be corrected to conserve 
mass (equation (94)). This is done such that the residuals of equation (95) 
remain unchanged. Therefore, the corrections in the second step are governed by 
the following discrete equations 

* - (V h *q - s) (96) 

-A h 5q + v h fip = o (97) 

where and are the discrete gradient and Laplacian operators on the given 
grid. 

To solve equations (96) and (97), one can introduce a discrete potential 
function. 

fiq = v h * (98) 


Therefore, 


V = 6 p = - (v h *q - s) (99) 

The discrete equation (99) is used to update the pressure and the velocity fields. 

In the first step of such calculations, the vorticity field is established. The 
potential correction in the second step does not alter the vorticity field, but it 
enforces the conservation of mass. 

These ideas are as old as the pressure correction methods for the solution of 
Navier Stokes equations of incompressible flows, introduced by Harlow and Welsh 
(48), Chorin (49), Temam (50), Patankar and Spalding (51), and others. Kim and 
Moin (52) used a similar decomposition for the continuous problem. Obviously, the 
discrete potential is preferred to avoid the difficulties associated with the 
boundary conditions for the intermediate variables. 

Returning to the transonic flow problem, one may consider the nonlinear Cauchy- 

Riemann problem 


(pu) x + ( P u) « o 


(100) 


330 



- w 


( 101 ) 


where p{ u, v) . 

The correction to the velocity field can be split into two parts, irrotational 
and rotational. These components can be calculated in terms of potential and 
stream functions. An important issue in this discussion is the treatment of 
shocks (and wakes). Fitting procedures are needed to calculate the entropy jump 
and the entropy gradient behind curved shocks. 

On the other hand, if the momentum equations are used to update the velocity 
field, and conservation of mass is enforced by a potential correction, it is not 
clear how to correct the pressure. Recently, there are some attempts to use the 
pressure correction methods for compressible flow calculations. The relation 
between the pressure correction and the velocity correction is, however, 
artificial. Still convergent results may be obtained for steady problems. 


VISCOUS FLOW SIMULATION 

In the previous sections, the velocity/vorticity formulation and the pressure 
correction methods for the solution of the Navier Stokes equations of 
incompressible viscous flows are discussed. Direct applications of vector 
potential methods are examined next. For example, in ref. (53), the following 
decomposition is used: 


u - + v, ■ v - v v, (I02) 

The continuity equation gives 

(p* x ) x + ( p ^y)y = 0 U 03 ) 

with the boundary condition 

*n ' *s/ P 004) 

where n and s are the normal and the tangential direction to the surface. The 
equation for 5 is 


<Wx + <Wy = - " 

(105) 

with the boundary condition 


v, ■ • 

(106) 


where w is obtained from the vorticity transport equation. It is assumed that the 
vorticity outside the viscous layer is negligible and hence a potential model is 
adequate. 

The boundary conditions (104) and (106) represent the coupling between a viscous 
and an inviscid problem. The perturbative stream function behaves like a 
displacement thickness which derives the potential calculations. On the other 
hand, the pressure distribution (<t> x ) derives the viscous calculations through the 
no slip boundary condition. No boundary-! ayer approximations have been made and 
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the formulation can be completely equivalent to the Navier Stokes equations. Only 
an inviscid grid is needed for the potential calculations, while the system is 
restricted to the viscous layer. Therefore, this formulation is a special form of 
domain decomposition techniques. 

Further approximations lead to some simplifications. For example, the density in 
the potential equation (103) can be calculated in terms of formula (V) and 
therefore the pressure, according to the isentropic relation is 

P - P?h M» 2 (107) 

On the other hand, the density in the system can be evaluated from the total 
enthalpy relation 

H ‘ + 2 + v2 > (108) 

In the transonic regime, H can be assumed constant even in the viscous layer, 
otherwise H can be evaluated from the energy equation. Unlike p., the density p 
has a boundary-layer- type profile like the velocity and the temperature. The 
pressure, however, can be assumed to be the same in the inviscid and viscous 
calculations. 

One way to implement this formulation is to use, in the viscous layer, the full 
*-a> system, with the surface boundary conditions, *«o and «o. The outer boundary 
conditions are w=o and ♦ * p<f> . Equation (108) is used for the density where the 
pressure is assumed knowrv. 

The next step is to construct 5. Equation (105) is solved with the surface 
boundary conditio^ given by (106). The boundary condition at the edge of the 
viscous layer ijs 5 = o and w is known from the first step. The output of this 
calculation is $ on n the surface. 

Finally, the potential equation is solved and the process is repeated until 
convergence. The problem for 5 in the second step may be solved coupled with the 
potential equation to avoid convergence difficulties. In this case, the *-w system 
provides u> to the 5-^ system, while the latter feeds back the pressure field and ^ 
at the edge of the viscous layer. s 

Extension to three-dimensional flows are possible using multiple stream 
functions. Alternatively, the viscous calculations can be based on the 
velocity/vorticity formulation. The role of the potential is simply to provide an 
approximation for the pressure field, which is needed to calculate the density in 
compressible flow cases. For incompressible flows, the potential can be dropped, 
since the pressure is eliminated from both the stream function and the 
velocity/vorticity formulations. The only advantage of keeping the potential is to 
restrict the domain of the viscous calculations. 

Recently, El Dabagbi (54) obtained finite-element solutions of Navier Stokes 
equations using a rotational correction to the potential field. He used the 
decomposition 


q = V0 + q, where v*q = o 


(109) 
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He solved the following potential equation 


(110) 


with the surface boundary condition <f> * o. The equations for q are very similar 
to incompressible Navier Stokes equations. In this formulation, the grid for 4 , 
problem must be fine enough to allow for the resolution of the nonhomogeneous term 
Vp* q. Similar ideas were introduced earlier by Dodge (55), Briley (56), Dwoyer 
(57) and others. 

Also, Ecer et al . (58) extended their formulation to treat viscous flows. Their 
approach is applicable for inviscid flow in the limit of high Reynolds numbers. 
Also, outside the viscous layer, the potential formulation is recovered. 

A composite velocity procedure for potential, Euler and Navier Stokes equations 
was introduced by Rubin et al . (59). The following velocity decomposition was 
proposed. 

u - (U + l)* x (111) 

v * * y (112) 

The multiplicative composite form of the axial velocity component consists of the 
potential component modified by the viscous velocity. The continuity and the 
tangential momentum equation are solved in a coupled manner to update <f> and U. 
The normal momentum equation in nonconservative form, is used to calculate the 
entropy. Numerical results for nonlifting airfoils indicate the differences 
between potential, Euler and Navier Stokes solutions. 

Helmholtz decomposition, vorticity generation and trail ing-edge condition for 
incompressible inviscid and viscous flows are discussed by Morino(60), (61). He 
derived also boundary integral equations for unsteady vjscous^ and inviscid flows. 
In this work, the solution of the Poisson's equations, « -Z is given by 

5 -fc < 113 > 

The rotational component of the velocity vector is then 

< 114 > 

The integral equation approach is also adopted by Kandil and Yates (62) to 
calculate inviscid vortical transonic flows, where both shocks and wakes are 
fitted. 

In this regard, the work of Wu (63), on numerical boundary conditions for 
viscous flow problems should be mentioned. Applications to multiple-body problems 
are given in ref. (64) . 
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CONCLUSIONS 


Unfortunately, no satisfactory results, based on vector potential methods, are 
available, so far, for real aerodynamic problems. For example, the simulation of 
transonic flow over a wing, with proper treatment of shocks and wakes, is still 
missing. At the same time, efforts on conservation laws gain momentum from recent 
mathematical developments in this field. Nevertheless, one can say, there is a 
great potential in the alternate formulations, particularly if practical methods 
for shock (and wake) fitting* are developed. 
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